function [fidSpectrum, fft_data] = GET_FIDSpectrum(obj, fid_signal, varargin)
%GET_FIDSPECTRUM Summary of this function goes here
%   Detailed explanation goes here
    p=inputParser;
    p.addParameter('func', @abs, @(x) isa(x, 'function_handle') );
    p.addParameter('isotope', 'both', @(x) any(validatestring(x, {'129', '131', 'both'})) );
    p.addParameter('cutOffFreq', 5.0, @(x) isscalar(x) && x > 0.0);
    p.addParameter('minProminence', 0.01, @(x) isscalar(x) && x > 0.0);
    p.addParameter('minDistance', 10.0, @(x) isscalar(x) && x > 0.0);
    p.addParameter('isPlot', false, @islogical);
    p.addParameter('signal_channel', obj.controller.DemodSettings.ByDemodFast.index, @(x) ismember(x, {'0', '1', '2', '3'}));
    p.addParameter('signal_component', 'x', @(x) ismember(x, {'x', 'y'}));
    p.addParameter('cut_factor', 10, @(x) isscalar(x) && x > 0);
    p.addParameter('nfft', 24, @(x) isscalar(x) && x>20);
    p.parse(varargin{:});

    f129est = obj.getXe129EstFreq();
    f131est = obj.getXe131EstFreq();
    fidData = fid_signal.getNodeData(p.Results.signal_channel, {p.Results.signal_component});
    fft_data = fftData(fidData, p.Results.nfft);
    [peak129, peak131] = locate_peak_raw(fft_data, f129est, f131est, p.Results.cut_factor);
    resonance129 = fitPeak(peak129);
    resonance131 = fitPeak(peak131);
    fidSpectrum=cell(1, 2);
    fidSpectrum{1}.resonance = resonance129;
    fidSpectrum{2}.resonance = resonance131;
    
    
%     switch p.Results.isotope
%         case 'both'
%             nPeak = 2;
%         case {'129', '131'}
%             nPeak = 1;
%     end
%     
%     fidData = fid_signal.getNodeData(p.Results.signal_channel, {p.Results.signal_component});
%     try
%         samplingRate=1/(fidData.time(2)-fidData.time(1));
%     catch
%         error('Cannot get sampling rate.');
%     end
%     
%     fidData.value=double([fidData.value;zeros(2^21-length(fidData.value),1)]);
%     fidData.value = fidData.value-mean(fidData.value);
%     points=length(fidData.value);
%     fft_data.frequency=samplingRate*(0:points/2)/points;
% 
%     amp = 2*p.Results.func( fft(fidData.value) )/points;
%     fft_data.amplitude=amp(1:points/2+1);
%     
%     normalizeData=fft_data.amplitude/max(fft_data.amplitude);
%     idx = find(fft_data.frequency > p.Results.cutOffFreq);    
%     [~,locs,w,~] = findpeaks(normalizeData(idx), fft_data.frequency(idx), ...
%                             'SortStr', 'descend', ...
%                             'Npeak', nPeak, ...
%                             'MinPeakDistance', p.Results.minDistance, ...
%                             'MinPeakProminence', p.Results.minProminence);
% 
%     fidSpectrum=cell(1, nPeak);
%     for k = 1:nPeak
%         fft_k=cut_spectrum(fft_data, locs(k), w(k), p.Results.cut_factor);
%         fidSpectrum{k} = getFIT(fft_k, 'isPlot', p.Results.isPlot);        
%     end

end

function res=fftData(signal, nPad)
    samplingRate=1/(signal.time(2)-signal.time(1));
    value=double([signal.value; zeros(2^nPad-length(signal.value),1)]);
    points=length(value);
    frequency=samplingRate*(0:points/2)/points;
    amp = 2*abs( fft(value) )/points;
    amplitude = amp(1:points/2+1);
    res.frequency = frequency(:);
    res.amplitude = amplitude;
end

function [res129, res131] = locate_peak_raw(fftdata, f129est, f131est, factor)
    lowFreqCutIdx = find(fftdata.frequency>5.0);
    freq = fftdata.frequency(lowFreqCutIdx); amp = fftdata.amplitude(lowFreqCutIdx);
    normAmp = amp/max(amp);
    
    [~,locs,w,~] = findpeaks(normAmp,  freq, ...
        'SortStr', 'descend', ...
        'Npeak', 2, ...
        'MinPeakDistance', 0.5*abs(f129est-f131est), ...
        'MinPeakProminence', 0.01);

    idx1 = find(fftdata.frequency < locs(1)+factor*w(1) );
    idx2 = find(fftdata.frequency > locs(1)-factor*w(1) );
    idx = intersect(idx1, idx2);
    res129.frequency = fftdata.frequency(idx);
    res129.amplitude = fftdata.amplitude(idx);
    res129.width = w(1);
    res129.w0 = locs(1);

    idx1 = find(fftdata.frequency < locs(2)+factor*w(2) );
    idx2 = find(fftdata.frequency > locs(2)-factor*w(2));
    idx = intersect(idx1, idx2);
    res131.frequency = fftdata.frequency(idx);
    res131.amplitude = fftdata.amplitude(idx);
    res131.width = w(2);
    res131.w0 = locs(2);

end

function [vMax, f0, width] = characterize_peak(peak)
    [vMax, idx] = max(peak.amplitude);
    f0=peak.frequency(idx);
    peakFun = @(x) interp1(peak.frequency,peak.amplitude,x) - vMax/sqrt(2.0);
    
    flist = min(peak.frequency):0.00001:max(peak.frequency);
    vlist = peakFun(flist);
    [delta, zIdx] = min(abs(vlist));
    
    if delta > vMax*1e-3
        error('zero point not found');
    end
    width = abs(flist(zIdx)-f0);
end

function resonance=fitPeak(peak)
    [vMax, f0, width] = characterize_peak(peak);    
    [xData, yData] = prepareCurveData( peak.frequency, peak.amplitude );

    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    
    
    ft=Constant.FittingModel('modelName', 'SqrtLorentz');
    opts.Lower =      [0.9*vMax*sqrt(pi)  0.0        0.99*width f0-0.1]; % a, b, gamma, w0
    opts.Upper =      [1.1*vMax*sqrt(pi)  0.0        1.01*width f0+0.1];
    opts.StartPoint = [vMax*sqrt(pi)      0.0        width      f0];
    opts.TolX = 1e-10;
    
    [res,gof] = fit( xData, yData, ft, opts );
    
    resonance.frequency = res.w0;
    resonance.gamma = res.gamma;
    resonance.hight =  res.a-res.b;
    resonance.data = peak;
    resonance.fit = res;
    resonance.gof = gof;
end


%% deprecated
function fidSpectrum = getFIT(fftData, varargin)
    p=inputParser;
    p.addParameter('isPlot', false, @islogical);
    p.parse(varargin{:});

    [xData, yData] = prepareCurveData( fftData.frequency, fftData.amplitude );
    [vMax, idx] =max(abs(yData));           
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';


    ft=Constant.FittingModel('modelName', 'SqrtLorentz');                    
    opts.Lower =      [0.0  0.0  0     min(fftData.frequency)]; % a, b, gamma, w0
    opts.Upper =      [Inf Inf 1.0   max(fftData.frequency)];
    opts.StartPoint = [vMax 0.0  0.010 fftData.frequency(idx)];

    [fidSpectrum.result, fidSpectrum.gof] = fit( xData, yData, ft, opts );

    fidSpectrum.frequency=fftData.frequency;
    fidSpectrum.amplitude=fftData.amplitude;
    
    fidSpectrum.resonance.frequency = fidSpectrum.result.w0;
    fidSpectrum.resonance.gamma = fidSpectrum.result.gamma;
    fidSpectrum.resonance.hight =  fidSpectrum.result.a-fidSpectrum.result.b;
    fidSpectrum.resonance.data = fftData;
    fidSpectrum.resonance.fit = fidSpectrum.result;
    fidSpectrum.resonance.gof = fidSpectrum.gof;
    
    if p.Results.isPlot
        fig = figure('Name', 'FID Spectrum'); ax = axes(fig);
        plot(fidSpectrum.result, fftData.frequency,fftData.amplitude);
        grid(ax, 'on');
        xlabel(ax, 'Frequency (Hz)');
        ylabel(ax, 'Normalized Amplitude');
        ylim(ax, [0, 1.1*vMax]);
    end
end

function res = cut_spectrum(spect, loc, width, factor)
    if nargin == 3
        factor = 20;
    end
    idx1 = find(spect.frequency < loc+factor*width);
    idx2 = find(spect.frequency > loc-factor*width);
    idx = intersect(idx1, idx2);
    res.frequency = spect.frequency(idx).';
    res.amplitude = spect.amplitude(idx);
end
